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ABSTRACT 

We introduce 21CMMC: a parallelized, Monte Carlo Markov Chain analysis tool, in¬ 
corporating the epoch of reionization (EoR) seminumerical simulation 21CMFAST. 
21CMMC estimates astrophysical parameter constraints from 21 cm EoR experiments, 
accommodating a variety of EoR models, as well as priors on model parameters and 
the reionization history. To illustrate its utility, we consider two different EoR sce¬ 
narios, one with a single population of galaxies (with a mass-independent ionizing 
efficiency) and a second, more general model with two different, feedback-regulated 
populations (each with mass-dependent ionizing efficiencies). As an example, com¬ 
bining three observations (z = 8, 9 and 10) of the 21 cm power spectrum with a 
conservative noise estimate and uniform model priors, we find that interferometers 
with specifications like the Low Frequency Array/Hydrogen Epoch of Reionization 
Array (HERA)/Square Kilometre Array 1 (SKAl) can constrain common reioniza¬ 
tion parameters: the ionizing efficiency (or similarly the escape fraction), the mean 
free path of ionizing photons and the log of the minimum virial temperature of star¬ 
forming haloes to within 45.3/22.0/16.7, 33.5/18.4/17.8 and 6.3/3.3/2.4 per cent, ~ ler 
fractional uncertainty, respectively. Instead, if we optimistically assume that we can 
perfectly characterize the EoR modelling uncertainties, we can improve on these con¬ 
straints by up to a factor of ~few. Similarly, the fractional uncertainty on the average 
neutral fraction can be constrained to within < 10 per cent for HERA and SKAl. By 
studying the resulting impact on astrophysical constraints, 21CMMC can be used to 
optimize (i) interferometer designs; (ii) foreground cleaning algorithms; (iii) observing 
strategies; (iv) alternative statistics characterizing the 21 cm signal; and (v) synergies 
with other observational programs. 

Key words: galaxies: high-redshift - intergalactic medium - cosmology: theory - 
dark ages, reionization, first stars - diffuse radiation - early Universe 


1 INTRODUCTION 


The epoch of reionization (EoR) describes the period of en¬ 
lightenment following the cosmic dark ages, when density 
perturbations grow into the first astrophysical sources (stars 
and galaxies). These sources produce ultraviolet (UV) ion¬ 
izing photons, which escape into the intergalactic medium 
(IGM) eventually reionizing the pervasive neutral hydrogen 
fog. 


The EoR is rich in astrophysical information, probing 
the formation and evolution of structure in the Universe, 
the nature of the first stars and galaxies and their impact 


on the IGM (see e.g. Barkana fc Loeb|2007 Loeb & Furlan 


|etto||2013| |Zaro ubi 2013). Currently the astrophysics of the 

EoR are poorly understood. However, a wave of upcoming 


* E-mail: bradley.greig@sns.it 


observations should produce a rapid advance in our under¬ 
standing. At the forefront of these will be the detection of 
the redshifted 21 cm spin-flip transition of neutral hydrogen 


(see e.g. Furlanetto et al. 2006b Morales & Wyithe 2010 
Pritchard & Loeb 2012) from several dedicated radio inter¬ 


ferometers. 

The first-generation 21 cm experiments such as the 


Yatawatta et al. 2013 

q the Murchison Wide Field Array 

(MWA; Tingay et al.| 

2013|Z 

and the Precision Array for 


Probing the Epoch of Reionization (PAPER; Parsons et al. 
2010jf] are just now coming online. These have enabled 


progress in understanding and characterizing the relevant 


1 http://www.lofar.org/ 

2 http://www.mwatelescope.org/ 

3 http://eor.berkeley.edu/ 
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2 B. Greig et al. 


systematics and foregrounds, and may even yield a detec¬ 
tion of the 21 cm power spectrum (PS) during the advanced 
stages of reionization (e.g. Chapman et al. 2012 Mesinger 
et al.|20T4 1 . 


In the near future, second- generation exper iments, the 
Square Kilometre Array (SKA; Mellema et al. 2013|p and 


the Hydrogen Epoch of Reionization Array (HERA; Beards- 
ley et al. 2015 will feature significant increases in collecting 


area and overall sensitivity. In addition to statistical mea¬ 
surements of the 21 cm PS, these second-generation experi¬ 
ments should provide the first tomographic maps of the 21 


cm signal from the EoR (Mellema et al. 2013s, deepening 


our understanding of reionization-era physics. 

However, we are faced with a fundamental question, 
what exactly can we learn from these observations? Quali¬ 
tatively, we possess several key insights. Reionization driven 
by relatively massive galaxies should result in larger, more 
uniform ionized regions (e.g. McQuinn et al. 2007| Iliev et al. 


2012) and the large-scale PS should peak at the midpoint of 
the EoR (e.g.|Lidz et al.|2008[), unless reionization is driven 


by X-rays (Mesinge r et al.|2013 l. Abundant, absorption sys¬ 
tems should result in an extended reionization (e.g. Ciardi 
|et al.||2006[) character ized by small, disjoint ionized regions 
(e.g. Alvarez &; Abel||2012 l. Although these studies provide 
valuable intuition, they do not quantify the constraints and 
degeneracies among astrophysical parameters. 

Therefore, it is imperative to develop an EoR. analy¬ 
sis tool to quantitatively assess what astrophysical informa¬ 
tion can be gleaned from the cosmic 21 cm signal. This will 
serve to guide current and future 21 cm experiments in opti¬ 
mizing their observing strategies and interferometer designs 
to maximize their EoR scientific returns. In recent years, 
several authors have explored the benefits of this approach 
by performing EoR parameter searches. However, these ap¬ 
proaches have been restricted to (i) sampling a finite, fixed 
grid of either simple analytic models (Choudhury & Fer- 
rara||2005 Barkana]|2009 1 or 3D seminumerical simulations 


(Mesinger et al. 2012 Zahn et al. 2012[ Mesinger et al. 2013); 


(ii) performing a Fisher Matrix analysis ( |Pober et al. 20141; 

or (iii) Monte Carlo Markov Chain (MCMC) sampling only 


the reionization history in an analytic framework (Harker 


et al.|20T2 Morandi fe Barkana|2012{|Patil et al.|2014l 

Here, we improve on these approaches by develop¬ 
ing a public 21 cm EoR analysis tool, named 21CMMC. 
21CMMC uses an optimized version of the well-tested, pub¬ 


lic simulation tool, 21CMFAST (Mesinger & Furlanetto 2007 
Mesinger et al.|[2011|), to generate 3D tomographic maps of 


the 21 cm brightness temperature, statistically comparing 
these to a (mock) observation. It quantifies the constraints 
and degeneracies of astrophysical parameters governing the 
EoR, within a full Bayesian MCMC framework. 

The remainder of this paper is organized as follows. In 
Section [2] we outline the development and implementation 
of 21CMMC. In Section [3j we highlight its performance at 
providing astrophysical constraints from a popular theoreti¬ 
cal model of the EoR containing a single population of ioniz¬ 
ing galaxies with a mass-independent ionizing efficiency. To 
further emphasize the strength and flexibility of 21CMMC 


4 https: //www.skatelescope.org 

5 http://reionization.org 


in Section [4j we then estimate the recovery of astrophysical 
parameters from a more generalized EoR model contain¬ 
ing two ionizing galaxy populations characterized by their 
mass-dependent ionizing efficiency. Finally, in Section [5] we 
summarize the performance of 21CMMC and finish with our 
closing remarks. Throughout this work, we adopt the stan¬ 
dard set of ACDM cosmological parameters: (I7 m , Ha, fib, 
n, a s , Ho) =(0.27, 0.73, 0.046, 0.96, 0.82, 70 kms" 1 Mpc -1 ), 
measured from the Wilkinson Microwave Anisotropy Probe 


1 Bennett et al. 2013 

1 which are consistent with the latest 

results from Planck 1 

Planck Collaboration XVIJ2014). 


2 21CMMC 


Prior to showcasing the performance of 21CMMcJ^] we first 
outline the various ingredients that contribute to the overall 
analysis tool. First, in Section we summarize the in¬ 
cluded EoR simulation code 21CMFAST. We then discuss 
the likelihood statistic in Section [2.2| Next, in Section [2.3| 
we briefly describe the specifics of the MCMC driver and in 
Section 2.4 we discuss the telescope sensitivities. Finally, in 
Section 2.5 we summarize the full 21CMMC pipeline. 


2.1 21CMFAST 


21CMMC employs a streamlined version of 21CMFAST_vl.l, 
a publicly available semi-numerical simulation code specifi¬ 
cally designed for enabling astrophysical parameter searches. 
It employs approximate but efficient methods for EoR 
physics, and is accurate compared to computationally ex¬ 
pensive radiative transfer (RT) simulations on scales rele¬ 
vant to 21 cm interferometry, ^ 1 Mpc (Zahn et alTpOll I. 

Specifically, 21CMFAST generates the IGM density, ve¬ 
locity, source and ionization fields by first creating a 3D real¬ 
ization of the linear density held within a cubic volume and 
then evolving the density held using the Zel’dovich approxi¬ 
mation (ZePdovich 1970) before smoothing on to a lower res¬ 
olution grid. Following this, 21CMFAST estimates the ion¬ 
ization held by comparing the time-integrated number of 
ionizing photons to the number of baryons within spherical 
regions of decreasing radius, R, following the excursion-set 


approach described in Furlanetto et al. (20041. A cell at co¬ 


ordinates ( x , z) within the simulation volume is then tagged 
as fully ionized if, 


Cfcoll (*^5 Z, R, -M m i n ) (5- 1, 


(1) 


where / co ii(a:, z, R, M m i n ) is the fraction of collapsed mat¬ 
ter within a spherical radius R residing within haloes larger 


than M min ( 

Press & Schechter 1974 

Bond et al. 1991 

Lacey 

& Cole 1993 Sheth & Tormen 

19991 and is an ionizing ef- 


hciency describing the conversion of mass into ionizing pho¬ 


tons (see Section 2.1.11. Partial ionization is additionally 
included for cells not fully ionized by setting the ionized 
fraction of a cell smoothed at the minimum smoothing scale, 

Rcell to f (.SC. Z. /V.cel], Miniri)- 


6 21CMMC will be made available at 
http://homepage.sns.it/mesinger/21CMMC.html Interested 

users can obtain a pre-release beta version by contacting the 
authors. 
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21CMMC: astrophysics from the 21 cm EoR signal 3 


It is common to characterize the EoR with just three pa¬ 
rameters: (i) the ionizing efficiency of high -2 galaxies; (ii) the 
mean free path of ionizing photons; (iii) the minimum virial 
temperature hosting star-forming galaxies. These parame¬ 
ters are somewhat empirical and in most cases are averaged 
over redshift and/or halo mass dependences. Nevertheless, 
they offer the flexibility to describe a wide variety of EoR 
signals, and can have a straightforward physical interpreta¬ 
tion. Below, we summarize each of these three parameters 
and their physical origins. 


2.1.1 Ionizing efficiency, 

The ionizing efficiency of high -2 galaxies (equation [T| can 
be expressed as 

<- 30 (tS) (m) (iSo 


1 + n T , 


( 2 ) 


where, f esc is the fraction of ionizing photons escaping into 
the IGM, /* is the fraction of galactic gas in stars, 1V 7 is the 
number of ionizing photons produced per baryon in stars 
and n Tec is the typical number of times a hydrogen atom 
recombines. While our reionization model only depends on 
the product of equation , we provide reasonable estimates 


for the individual factors. The choice of N-, 


4000 is ex¬ 


pected from Population II stars (e.g. Barkana fe Loeb|2005 1, 
while /+ and / eS c are extremely uncertain in high -2 galaxies 


(e.g. Gnedin et al.|2008 Wise fe Cen||2009 Ferrara & Loeb 


20131, though our choices are consistent with observations 


of galaxy luminosity functions (e.g. Robertson et al. 2013 


Cooke et al. 2014 Dayal et al. 2014 1 . Finally, n 


consistent with the models of (Sobacchi & Mesinger 20141 


1 is 


which result in a ‘photon-starved’ end to reionization, con¬ 
sistent with emissivity estimates from the Lya forest (e.g. 
|Bolton fe Haehnelt||2007| |McQuinn et al.||2011| |. We allow 
£ to vary within the range C, £ [5,100]. Since, (j is empiri¬ 
cally defined unlike / eac which can be observationally con¬ 
strained at lower redshift, we will often provide both when 
recovering our parameter constraints. For our range of £, 
this corresponds to / esc £ [0.05,1] using our fiducial values 
adopted in equation ( 2j. In this work, we will first consider a 
constant (Section[3l before considering a more generalized 
halo mass-dependent model in Section [4] 


2.1.2 Mean free path of ionizing photons within ionized 
regions, R m t p 

The propagation of ionizing photons through the IGM 
strongly depends on the abundances and properties of ab¬ 
sorption systems (Lyman limit as well as more diffuse sys¬ 
tems), which are below the resolution limits of EoR simula¬ 
tions. These systems act as photon sinks, roughly dictating 
the maximum scales to which H ii bubbles can grow around 
ionizing galaxies. In EoR modelling this effect is usually 
parametrized with an implicit maximum horizon for ioniz¬ 
ing photons (corresponding to the maximum filtering scale in 
excursion-set EoR models). Physically, this maximum hori¬ 
zon is taken to correspond to the mean free path of ionizing 


photons within ionized regions, i?. m f p (e.g. O’Meara et al. 


2007| [Prochaska fe Wolfe|2009HSongaila fe Cowie|2010||Mc 


Quinn et al. 2011). Motivated by recent sub-grid recombi¬ 
nation models (Sobacchi & Mesinger 20141, here we explore 
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the range R m f p £ [5, 20] cMptQ While this range is relatively 
narrow, we point out that this parameter only becomes im¬ 
portant once the size of the ionized regions become larger 
than J? m f p (e.g. McQuinn et al.||2007 Alvarez fe Abel|2012 


Mesinger et al.|2012 1. Expanding the allowed range of R m f P 

we consider in this work would therefore not greatly modify 
our conclusions. 


2.1.3 Minimum virial temperature of star-forming haloes, 

rji Feed 
vir 

Throughout, we choose to define the minimum threshold 
for a halo hosting a star-forming galaxy to be in terms of 
its virial temperature, which regulates processes important 
for star formation: gas accretion, cooling and retainment of 
supernovae outflows. The virial temperature is related to the 
halo mass via, (e.g. Barkana fe Loeb]|2001 1 

(&)-'*'" A '- /a 


TV i 


1.98 x 10 4 K 


Qj 

3/2 


A 


187r 2 y 

1 + 2 
10 


- 3/2 


Mrr 


0’ 


(3) 


where fi is the mean molecular weight, = £2m(l + 
z) 3 /[O m (l + z) 3 + n A ], and A c = 18 tt 2 + 82 d — 39 d 2 where 
d = — 1. Typically, T v - ir « 10 4 K has been adopted in 

the literature as the minimum temperature when efficient 
atomic cooling occurs, corresponding to a DM halo mass of 
~ 10 8 Mq at z ~ 10. In principle, this could be lower as the 
first star s are likely hosted within haloes with Mhaio ~ 10 6 - 
10 7 Mq ( Haiman et al.|l996 


Abel et al 


12002 ] 

these h 


Bromm et al. 


2002). However, star formation within these haloes is likely 
inefficient (a few stars per halo) and can be quickly (z > 20) 
suppressed by Lyman-Werner or other feedback processes 


well before the EoR i 

Haiman et al.|2000 

Ricotti et al. 2001 

Haiman & Bryan 2006 Mesinger et al. 2006 Holzbauer & 

Furlanetto|2012 

Fialkov et al. 2013 

)• 

Here, we allow for a critical 

temperature threshold, 


which can even be larger than 10 4 K. Below T^ e 


7 The ionization structure of these subgrid models of recombi¬ 
nations inside unresolved systems does not in fact directly trans¬ 
late to a constant, uniform value of R m f p . These authors find 
that it takes a long time for cosmic Hu regions to approach 
the Stromgren limit, with their expansion only slowing more and 
more as reionization progresses. In other words, the mean free 
path (which is only defined in terms of the instantaneous re¬ 
combination rate) is almost always larger than the actual size 
of the Hu regions (which depend on the cumulative number of 
recombinations), even at the late stages of reionization. The re¬ 
sulting ionization fields have a dearth of ionization structure on 
scales smaller than the instantaneous mean free path. Moreover, 
the spatial correlations between the sources and sinks of ioniz¬ 
ing photons results in sizeable spatial and temporal variations in 
the mean free path. Nevertheless, an effective, uniform maximum 
horizon for ionizing photons of ~ 10 Mpc can reproduce the 21 
cm power spectra of the |Sobacchi Sz Mes inger (2014j simulations 
to within tens of percent on relevant scales, at the epochs studied 
in that work 0.25 < rm < 0.75 (Sobacchi, private communica¬ 
tion; Mesinger Sz Greig, in preparation). For historical context, 
we use the variable ‘-Rmfp’ to denote this effective horizon set by 
subgrid recombinations, but caution that the actual mean free 
path is larger than this parameter, as discussed above. 























































































4 B. Greig et al. 


(but above 10 4 K), haloes are still sufficiently large to pro¬ 
duce low-mass galaxies; however, their potential wells are 
not deep enough to retain supernovae-driven outflows that 
quench star formation and their corresponding contribu¬ 


tion to reionization (e.g. Springel & Hernquist 2003). In 
this work, we allow T^® ed to vary within the range T *;® ed £ 
[10 4 , 2 x 10 5 ] K, with the lower limit set by the atomic cool¬ 
ing threshold and the upper limit consistent with observed 
high-z Lyman break galaxies (see for example Fig. [lj. 


2.2 Likelihood statistic 

The offset of the 21 cm brightness temperature relative to 
the CMB temperature, T 7 (e.g. Furlanetto et al.|2006b I, can 
be written as, 


ST b (u) = 


Ts T "' (l-e^) 
1 + z v ' 

27XHl(l + <5nl) 


H 


1 + 2 0.15 

10 Q m h 2 


dur/d r + H 
1/2 / Q. h h i 


l-£ 

T s 


0.023 


mK, 


(4) 


where Ts is the gas spin temperature, r„ 0 is the optical depth 
at the 21 cm frequency, t'o, 5 n i(x,z) is the evolved (Eular- 
ian) overdensity, H(z) is the Hubble parameter, du r /dr is 
the gradient of the line of sight component of the velocity 
and all quantities are evaluated at redshift z = vq/v — 1. 
Throughout this work, we include the effects of peculiar ve¬ 
locities and assume Ts T 7 (likely appropriate for the ad¬ 
vanced stages of reionization, corresponding to an average 
neutral fraction of ihi < 0.8, e.g. Mesinger et al. 20131 


as the explicit computation of the spin temperature within 
21CMFAST is computationally expensive; however, in future 
we will remove this assumption. 

Although 21CMFAST produces 3D maps of the 
21 cm brightness temperature, we require a good¬ 
ness of fit statistic in order to quantitatively com¬ 
pare an EoR realization with the mock observation. 
We choose the common spherically averaged 21 cm 
PS, A h{k,z) = k 3 /(2n 2 V)5T 2 (z)(\8 2 i(k,z)\ 2 } k where 
82 i(x,z ) = STs(x, z)/5Tb(z) — 1, to be this statistic. It is 
important to note that while we restrict our default analy¬ 
sis to A 21 (k,z), in principle any statistical measure of the 
cosmic 21 cm signal could be used in our framework. 


2.3 The MCMC driver 

An MCMC based EoR tool is a computationally challeng¬ 
ing prospect since each step in the computation chain re¬ 
quires a 21CMFAST realization. For reference, computing a 
128 3 21 cm realization at a single redshift takes ~ 3s us¬ 
ing a single core. This is efficient, however, it is a case of 
diminishing returns if we attempt to improve the compu¬ 
tation time using the inbuilt multithreading of 21CMFAST. 
Therefore, constructing an MCMC analysis tool which sam¬ 
ples 21CMFAST utilizing existing serial MCMC algorithms 
such as Metropolis-Hastings (Metropo lis et al7||1953 Hast¬ 
ings! 1970| ), which require multiple long computation chains, 
would be too computationally inefficient. 

In recent years, several new, more efficient and parallel 


MCMC algorithms have been developed. Once such algo¬ 
rithm is the affine invariant ensemble sampler developed by 


Goodman & Weare (20101. This algorithm has since been 


modified and improved before being released as the publicly 
available python module emce e by|Foreman-Mackey et alT| 
(20131. Recently, Akeret et al. (20131 discussed the virtues 


of adopting these computational algorithms for applications 
of cosmological parameter sampling. These authors released 
an easy to implement, massively parallel, publicly available 
python module COSMOHAMMER which has already been ex- 
tensibly adopted for a wide variety of cosmological applica¬ 
tion. We therefore choose to build 21CMMC using a modified 
version of COSMOHAMMER which performs our MCMC reion¬ 
ization parameter sampling by calling 21CMFAST at each 
step in the computation chain to compute the log-likelihood. 


2.4 Telescope noise profiles 

To illustrate the performance of 21CMMC for current and 
future 21 cm experiments, we compute the telescope sensi- 


tivities using the python module 21CMSENStpj |Pober et al. 
|2013|[2014| ). In this section, we outline only the specifics and 
assumptions we make in producing our realistic telescope 
noise profiles, summarizing the key telescope parameters in 
Table [TJ and defer the reader to the more detailed discus 


sions within Parsons et al. (2012a I and Pober et al. (2013 
2014). 


The thermal noise PS is computed at each un-cell ac¬ 
cording to the following (e.g. [Morales 2005. McQuinn et al. | 
2006]|Pober et al.|2014| ), 


(5) 


A 

where X 2 Y is a cosmological conversion factor between ob¬ 
serving bandwidth, frequency and comoving distance units, 
$Y is a beam-dependent factor derived in Parsons et al. 
(20141, t is the total time spent by all baselines within a 
particular k mode and T sya is the system temperature, the 
sum of the receiver temperature, T rec , and the sky tempera¬ 
ture T a ky. For all telescope configurations considered in this 
work, we assume drift scanning with a total synthesis time 
of 6 h per night. While HERA can only operate in drift 
scan mode, LOFAR and SKA can perform tracked scans. 
In the future, we will investigate various observing strate¬ 
gies; however, drift scanning suffices here for our like-to-like 
comparisons. 

Of great concern for 21 cm experiments is the estima¬ 
tion and removal of the bright foreground emission. Most no¬ 
table of these are the chromatic effects due to how the inter¬ 
ferometric array’s uv coverage depends on frequency. It has 
been shown that these effects do not lead to mode-mixing 
across all frequencies, but instead localize the majority of 
the foreground noise into a foreground ‘wedge’ within cylin- 


drical 2D fc-space ( 

Datta et al. 

2010 

Morales et al. 

2012 

Parsons et al. 2012b Trott et al. 2012 

Vedantham et al. 

2012 

Thyagarajan et al. 2013 

mu et al. 

2014a b 

1 . This re- 


sults in a relatively pristine observing window where the cos¬ 
mological 21 cm signal is dominated only by thermal noise. 


https: / / github. com/jpober/21cmSense 
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21CMMC: astrophysics from the 21 cm EoR signal 5 


Parameter 

LOFAR 

HERA 

SKA 

Telescope antennas 

48 

331 

866 

Diameter (m) 

30.8 

14 

35 

Collecting area (m 2 ) 

35 762 

50 953 

833190 

Tree (K) 

140 

100 

0-lTgky + 40 

Bandwidth (MHz) 

8 

8 

8 

Integration time (hrs) 

1000 

1000 

1000 


Table 1. Summary of telescope parameters we use to compute 
sensitivity profiles (see the text for further details). 


However, it is still uncertain where to define the transition 
from the foreground ‘wedge’ to the EoR observing window. 

Moreover, the summation over redundant baselines 
within an array configuration can reduce uncertainties in 
the thermal noise estimates (Parsons et al. 2012ai. However, 


the gain depends on whether they are perfectly or partially 


coherent (Hazelton et al. 

2013]). 

In Pober et al. 

(2014 

, these authors consider numerous 

foreground removal strategies, including a variety of treat- 


ments for the summation over redundant baselines and how 
to define the transition from the foreground ‘wedge’. We 
defer the reader to their work for more detailed discussions, 
and highlight that throughout this work, we adopt what they 
refer to as their ‘moderate’ foreground removal. For this, we 
assume a fixed location of the foreground ‘wed ge’ to extend 
Afcii = 0.1 h Mpc -1 beyond the horizon limit (Pober et al. 


2014) and compute our estimates of the thermal noise PS 


sampling only fc-modes that fall within the EoR observing 
window. This strategy further assumes the coherent addition 
of all redundant baselines. 

We include another level of conservatism to our errors, 
applying a strict fc-mode cut to our 21 cm PS. We sample 
only modes above |fc| = 0.15 Mpc -1 . Following the moder¬ 
ate foreground removal strategy of Pober et al. (20141, it 


effectively turns out that there is limited information avail¬ 
able below this fc-mode; however, we still apply this strict 
cut on top of the estimated total noise PS. Modes with |fc| ^ 
0.15 Mpc -1 can still be affected by the foreground ‘wedge’ 
due to its spherically asymmetric structure, hence we apply 
both cuts. 

Since equation (|5j) is computed for each individual k- 
mode, the sample variance of the cosmological PS can be 
easily included to produce an estimate of the total noise 
by performing an inverse-weighted summation over all the 
individual modes (jPober et al.|2013|), 

- 1/2 

2 ,,, /v- 1 ' (6) 


SA^ +S (k) — ( y 


(Al,(k) + AUk)) 2 


where <SA|. +s (fc) is the total uncertainty from thermal noise 
and sample variance in a given fc-mode and A^i (k) is the 
cosmological 21 cm PS. Here, we assume Gaussian errors for 
the cosmic-variance term, which is a good approximation on 
large scales. 

Within this work, we restrict our analysis to three tele¬ 
scope designs (LOFAR, HERA and the SKA). In jPober et al. 
(2014), these authors qualitatively compare the merits of a 


variety of telescope designs through a measurement of the 
overall sensitivity of the 21 cm PS. To perform this direct 
comparison, these authors set a variety of design specific 
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parameters to be equivalent. Instead, we perform our com¬ 
parison using the exact design specifics of each experiment 
as outlined below. Note, the difference is marginal, however 
using system specific parameters we better mimic the ex¬ 
pected overall sensitivities of the modelled telescope. For 
all, we model T s k y using th e frequency dependent scaling 


r sky = 60 ( 300 m H z ) 2 ' 5j K ( Thompson et al.|2007 |. 


1. LOFAR: we model LOFAR using the antennas po¬ 
sitions listed in |van Haarlem et al.| ([2013 ). Following the 
approach of Pober et al. (20141, we focus solely on the 


Netherlands core for the purposes of detecting the 21 cm PS 
assuming that all 48 high-band antennas (HBA) stations 
can be correlated separately maximizing the total number 
of redundant baselines. We assume T lec = 140 K consistent 
with I Jensen et al. ( [201 3). 

2. HERA: we follow the design specifics outlined in 


Beardsley et al. (2015), where they propose a final design 


including 331 antennas. Each antennas is 14m in diameter 
closely packed into a hexagonal configuration to maximize 


the total number of redundant baselines (Parsons et al. 


2012a). While the total observing area for HERA is an order 


of magnitude lower than the SKA, its design is optimized 
specifically for P S measure ments allowing for comparable 
sensitivity (Pober et al. [2014|^ However, our 331 antennas 
telescope is smaller than the larger 547 antennas instrument 
proposed by these authors, therefore, our design performs 
more as an intermediate instrument. For HERA, we model 
the total system temperature as T sya = 100 + T sky K. 

3. SKA: we mimic the current design brief for the SKA- 
low Phase 1 instrument outlined in the SKA System Baseline 
Design document! 19 ] Specifically, SKA-low Phase 1 includes 
a total of 911 35 m antennas stations, with 866 stations nor¬ 
mally distributed in radius within a core extending out to a 
radius of ~2 km and a further 45 stations arranged into three 
spiral arms each with 15 stations extending to ~ 50 km. 
For our SKA design we include only the core stations, as 
the spiral arms add little to the PS sensitivity. As in [Pober] 
et al. (2014), we model the central core of SKA as a hexagon 


to maximize baseline redundancy, however our core extends 
out to a radius of ~ 300m which includes a total of 218 
telescope antennas. We then randomly place the telescope 
antennas normalising to ensure we contain 433 and 650 tele¬ 
scopes within a radii of 600 and 1000 m and randomly place 
the remaining 216 antennas at a radii larger than 1000 m. 
The total SKA system temperature is modelled as outlined 
in the SKA System Baseline Design, T sys = l.lT sky -+- 40 K. 


9 These authors found a factor of 2.5 difference between their 
SKA and HERA designs for their estimated total sensitivity on 
a measurement of the 21 cm PS at z = 9.5. However, this factor 
was a result of a numerical error in their calculation of their SKA 
sensitivities. Accounting for this, their HERA design has in fact 
only a marginally reduced sensitivity relative to the SKA for their 
pessimistic and moderate foreground scenarios (Pober, private 
communication) 

10 http://www.skatelescope.org/wp- 
content/uploads/2012/07/SKA-TEL-SKO-DD-001- 
1 _BaselineDesignl .pdf 
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2.5 The 21CMFAST pipeline 

Eventually 21CMMC will be performed on actual observa¬ 
tions. However for this proof-of-concept, we first create mock 
observations of the cosmic EoR signal and generate sensitivi¬ 
ties for various instrument^] Using 21CMFAST, we perform 
a high-resolution, large-box realization of our fiducial EoR 
model, simulated within a 500 3 cMpc 3 box with an initial 
grid of 1536 s voxels smoothed on to a 256 s grid. We then 
calculate the 21 cm PS from this simulation at several dif¬ 
ferent redshifts and compute the total noise contribution 
from thermal and sample variance using 21CMSENSE. In ad¬ 
dition to the thermal noise and sample variance, we add 
an additional uncertainty due to EoR modelling, account¬ 
ing for, e.g. differences in the implementation of RT (e.g. 
Iliev et al. 2006I and errors in semi-numerical approxima¬ 
tions (e.g. Zahn et al. 2011). In other words, even if we have 
a well behaved Universe characterized exactly by our 3 or 5 
parameter models, the simulated PS will be different from 
the true one due to the numerical implementation. We as¬ 
sume for this a constant multiplicative error of 25 per cent, 
guided by the observed 10-30 per cent difference at the rel¬ 
evant scales between 21CMFAST and the cosmological RT 


simulations in Zahn et al. (2011). For each mock observation 


of the 21 cm PS, we add this 25 per cent error in quadrature 
with the noise contribution from equation §. In principle 
this would be a correlated error, however, at this stage a 
constant multiplicative factor is sufficient. Below, we include 
some results excluding this error, showing the maximum in¬ 
formation which can be extracted from the signal assuming 
a perfect characterization of modelling uncertainties. 

To sample the reionization parameter space in 
21CMMC, we perform a smaller, lower resolution EoR simu¬ 
lation. Using a different set of initial conditions to the mock 
observations, we simulate a 250 s cMpc 3 volume with an ini¬ 
tial grid of 768 3 voxels smoothed down on to a 128 3 grid. 
Our two simulations are both sampled at the same physi¬ 
cal resolution (~ 2 Mpc), and resampled by the same factor 
of 6 from the initial density grid to ensure convergence be¬ 
tween the recovered 21 cm PS from the mock observation 
and the PS from the smaller 128 3 simulations. To account 
for the impact of shot noise on the smallest spatial modes 
(largest k) within these low-resolution simulations, we apply 
an additional cut at k ^ 1.0 Mpc -1 . 

At each step in the MCMC algorithm we compute the 
21 cm PS corresponding to that point in our astrophysical 
parameter space using the 3D realization from 21CMFAST. 
We compute its likelihood using the y 2 statistic between 
k = 0.15 and 1.0 Mpc -1 to assess whether to accept or 
reject the new proposed set of EoR parameters. 

Finally, we briefly discuss the improvements in our 
method to sample the astrophysical EoR parameter space, 
relative to previous works. These previous approaches have 
been restricted to either sampling a fixed, course grid of 


reionization models (e.g. Barkana 2009 Mcsinger et al. 2012 
Zahn et al.||2012i or the fundamental assumption of Gaus- 


11 At present, 21CMMC computes only co-evolving data cubes of 
the EoR.. In future, we intend to improve the simulation pipeline 
to mimic a more realistic observational pipeline, taking into ac¬ 
count several effects such as finite light-cones, foreground cleaning 
etc. 


sian errors in Fisher Matrix applications (Pober et al. 2014 


although new methods exist to overcome this assumption, 
e.g. Joachimi fe Taylor|2011 Sellentin et al.|2014 |. Instead, 
with a Bayesian MCMC framework, we are able to more 
efficiently sample the EoR parameter space and by con¬ 
struction obtain direct estimates of the individual proba¬ 
bility distribution functions (PDFs) of the EoR parameters, 
removing the a priori assumptions on the intrinsic probabil¬ 
ity distribution (i.e. Gaussian). Moreover, by sampling over 
10 5 realizations within 21CMMC this allows us to sample 
many features within the EoR parameter space not visible 
in course-grid models. Additionally, this allows for the bet¬ 
ter characterization of parameter degeneracies, which is a 
problem for Fisher Matrix analyses. 


3 SINGLE IONIZING GALAXY POPULATION 


To illustrate the use of 21CMMC, we first adopt the popular 
three-parameter EoR model, discussed above. This model is 
characterized by a single population of efficient star-forming 
galaxies hosted by haloes above our defined critical temper¬ 
ature . The ionizing efficiency of this model is a step 
function, 


C(Tvir) = 


Co 

0 


1 • > T 

vir 

1 • <r T 

vir \ \ 


Feed 

vir 

iFeed 

vir 


(7) 


In this model, a constant fraction of the host halo mass 
goes into the production of ionizing photons. To illustrate 
the performance of 21CMMC, we choose £o = 30, R m f p = 
15 Mpc and T^-® ed = 3 x 10 4 K for our mock observation. 
Our chosen values are roughly consistent with observational 
data, though they merely serve for illustrative purposes. 


3.1 Building physical intuition 

3.1.1 UV luminosity function 


It is first constructive to illustrate how this EoR model 
would appear with respect to current observations. In Fig. [I] 
we provide a schematic picture of our single ionizing galaxy 
population compared to the observe d non - ionizing (~ 15 00 
A) UV luminosity function at z = 8 of Bouwens et al. (20141. 
To obtain the UV luminosity of our galaxy population, we 
use the Schechter function parameter values described in 


Kuhlen & Faucher-Giguere (2012) and abundance matching 


to estimate the host halo mass as a function of the UV mag¬ 
nitude, M(Muv). Fitting this expression as a power-law, we 
convert our T^ ed into a UV magnitude which produces a 
sharp drop in the UV luminosity function at Muv = —12.6 


(M 


10 s 


' Mq). By construction, the non-ionizing UV 


luminosity function matches the observational data and we 
extrapolate down into the faint-end until , below which 
no ionizing photons are produced/escape the galaxy. In our 
illustrative model, this drop, tracing fundamental galaxy for¬ 
mation physics, is well beyond the sensitivity limits even for 
JWST (e.g. Salvaterra et al.[2 0111, showcasing the power of 
our approach. 

The UV luminosity function with an appropriate dust 
model can constrain the star formation rate, which (approx¬ 
imately) collapses the uncertainty in the ionizing emissiv- 
ity to just that in the escape fraction, f esc . In other words, 
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logio(Total Halo Mass, M [Mq]) 

13.9 12.2 11.0 10.3 9.7 9.2 8.6 8.0 7.3 



Muv 


Figure 1. A schematic picture of the z = 8 non-ionizing UV 
luminosity function (blue curve) from our fiducial single popula¬ 
tion reionization model (used for the mock observation) char¬ 
acterized by a threshold temperature of TT® ed = 3 X 10 4 K 
(M ~ 10 8,78 Mq at jz = 8). Haloes above T^? ed host typical 
star-forming galaxies contributing to the reionization of the IGM 
(green shaded region). Below TT®, haloes are small enough to 
be affected by internal feedback processes (blue shaded region). 
At T v ir < 10 4 K (M rsj io 8 06 Mq at z — 8), gas can no longer 
efficiently cool to accrete on to haloes. Note that the model curve 
at T vir > T£f d is obtained from a power-law fit to the halo 
mass as a function of UV magnitude recovered from abundance 
matching, which we can use to imply that for our fiducial choice 
of a constant ionizing efficiency we roughly obtain the scaling, 
f esc oc M 0 - 32 (see the text for further discussions). 


the total ionizing luminosity, Li on is roughly proportional 
to the product of the non-ionzing UV luminosity (£uv) 
and /esc: L lon (M) oc f eac (M)Luv(M). Our fiducial choice 
of a constant £ implies Li on oc M, from which we can 
obtain an expression for Lw(M) using our power-law ex¬ 
pression for M(Muv) to recover L\jy(M) oc M°' 68 for the 
faint end of the luminosity function. Therefore, our illustra¬ 
tive, single population model at z = 8 roughly translates to 

/esc OC M°- 32 . 

3.1.2 How does the 21 cm PS depend on the EoR 
parameters? 

In Fig. [2] we highlight how each individual EoR parameter 
affects the 21 cm PS, and also provide a visual representation 
of the total noise profiles for each telescope. In all panels, 
we show the mock z = 9 observation for our fiducial EoR 
parameters (solid black curve) which results in ihi = 0.71. 
In order to facilitate the discussion, we decompose the 21 cm 
PS to first order (neglecting the contribution from redshift- 
space distortions and spin temperature fluctuations), 

A|r(fc) oc XniAss(k) - 2 x H i (1 - xhi) A] x (k) 

+(1-xhi) 2 A _ 2 x (fc), (8) 

where A 2 5 , A| x and A xx are the matter density, density- 
ionization and ionization PS respectively and fluctuations 
in the ionized fraction are defined as <5 X = (aimi — 

(*hii))/(*hii)). Although |Lidz et aTT| ( |2007| ) showed that this 
first-order approximation provides a rather poor description 


for the 21 cm PS, for our qualitative discussion below equa¬ 
tion Q sufficiently highlights the major components to the 
total power. 

Before reionization gets well underway (xhi > 0.9) A| x 
closely follows the matter density PS, A 2 5 . Initially, cos¬ 
mic H ii bubbles grow around biased density peaks hosting 
galaxies. This causes the large-scale power to drop due to the 
increased (negative) contribution from the cross-correlation 
term, A 2 X . As time evolves (xhi decreases), the Hu bubbles 
become more numerous (with the increasing population of 
ionizing sources) and begin to overlap, increasing the con¬ 
tribution from the ionization PS, A xx . Eventually, this term 
begins to dominate, resulting in an overall increase in the 
total power in A|i across all scales. At xhi = 0.5, the large- 
scale power peaks. As the Hu bubbles grow and overlap, the 
peak power in A 2 i shifts to larger physical scales (smaller 
k) as A xx is sensitive to the physical size of the H ii regions. 
At the same time, the power on small scales begins to drop 
as the number of small, isolated Hu regions decreases. 

This generic PS evolution during the EoR is evident in 
the first three panels of Fig. [2] (see also, e.g. McQuinn et al. 


|2007||Lldz et al.|2008||Friedrich et al.|2011[ ). A larger value of 
£o corresponds to a more advanced EoR at z = 9, as galaxies 
are more efficient at producing ionizing photons. Similarly, 
a smaller value of T^ ed results in a more advanced EoR, as 
smaller, more abundant galaxies can efficiently form stars. 
The progress of the EoR is less sensitive to variation in R m f p 
(second panel); instead this maximum horizon for the ion¬ 
izing photons produces a prominent ‘knee’ in the PS, sup¬ 
pressing coherent ionization structure on large scales (e.g. 
|Alvarez fe Abel|2012|[Mesinger et al.|2012[ ). 

The PS is not however entirely determined by xhi. In 
the final panel of Fig. [2] we illustrate the impact of our EoR 
parameters at a fixed IGM neutral fraction of aim = 0.71 
(i.e. all at the equivalent stage of reionization). For exam¬ 
ple, if the EoR is driven by more biased, larger mass haloes 
(larger Tfs ed ), their relative dearth must be compensated for 
by a higher ionizing efficiency. This results in an EoR mor¬ 
phology with lar ger, more isolated cosmic Hu patches (e.g. 
McQuinn et al.| [2007). A xx is more peaked on large scales, 


with a decrease in A^ x . As a result, the 21 cm PS is higher, 
with a more prominent large-scale ‘bump’ tracing the typi¬ 
cal H ii region size. Also explicitly evident in the last panel 
is the prominent ‘knee’ feature, imprinted by small values 
of -Rmfp. Here, for i? m f p = 3 Mpc (blue curve), the ioniza¬ 
tion PS peaks on much smaller scales, with a sharp drop in 
large-scale power compared to the other curves computed 
with R m f p = 15 Mpc. Therefore, the 21 cm PS contains 
information on both the EoR epoch (i.e. aim), as well as in¬ 
dependent (i.e. at fixed aim) information on the properties 
of galaxies and the IGM. 

In all panels, we provide the estimated la errors on 
the 21 cm PS for each of the considered 21 cm experi¬ 
ments, including the additional 25 per cent modelling un¬ 
certainty. It is immediately obvious that LOFAR will have 
difficulty discriminating among various models, having to 
rely on the largest scales close to the foreground-dominated 
regime. HERA and the SKA on the other hand should be 
able to recover meaningful constraints on EoR parameters, 
with the large-scales limited by the additional 25 per cent 
modelling uncertainty we include. On smaller spatial scales, 
SKA performs considerably better than HERA. 
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k (Mpc *) k (Mpc x ) 

Figure 2. The 21 cm PS at z = 9 for various values of our three reionization parameters, Qj. -R m f p and Tf® ed . In all panels, the thick 
black curve corresponds to our fiducial reionization model, with (£o, -R m f p , Tf eed ) = (30, 15 Mpc, 3 X 10 4 K) and a neutral fraction of 
Shi = 0.71. Light to dark shaded regions correspond to the 1 rj errors of the total noise PS (see Sections |2.4| and 1 2.5 ( for the three 21 
cm experiments we consider in this analysis (LOFAR, HERA and SKA), including a 25 per cent modelling uncertainty. Dashed vertical 
lines at k = 0.15 and k = 1.0 Mpc -1 demarcate the physical scales on which we choose to fit the 21 cm PS, while the hatched region 
corresponds to the foreground dominated region. We show the impact of the ionizing efficiency (y (top left), the maximum horizon for 
ionizing photons, f? m f p (top right) and the minimum virial temperature of star-forming galaxies, Tf® ed (bottom left). Bottom right: the 
21 cm PS corresponding to different EoR parameters, but at the same EoR epoch (ini) as the fiducial simulation. 


3.2 Single epoch observation of the 21 cm PS 


21 cm experiments have coverage over a large bandwidth 
(e.g. 50-350 MHz for the SKA allowing for observations to 
2 < 28). However, foregrounds and high data rates can limit 
the coverage to narrower instantaneous bandwidths. Hence, 
we begin by analysing an observation at a single redshift (as¬ 
suming our fiducial bandwidth of 8 MHz), before moving on 
to a broader bandpass observation. We restrict our analysis 
to only second-generation 21 cm experiments, HERA and 
the SKA, since with a single bandwidth the first-generation 
instruments will only achieve a marginal detection at best 
(e.g. Mesinger et al.|2014 Pober et al.|2014 see also Fig. pjj). 

In Fig. [3] we compare the outputs of 21CMMC for both 
HERA (red curve) and SKA (blue curve) for an assumed 
single z — 9 observation of our fiducial 21 cm PS with a 
total integration time of 1000 h. In this figure, and for the 
remainder of this work we assume uniform priors on all re¬ 


covered parameters within their allowed ranges outlined in 
Section [2.11 The dashed vertical lines denote the mock ob¬ 
servation values. Across the diagonal panels of this figure 
we provide the ID marginalized PDFs for each of our EoR 
parameters. In the top-right panel, we provide the marginal¬ 
ized ID PDFs of the IGM neutral fraction. Additionally, we 
choose to renormalize all ID PDFs to have the peak proba¬ 
bility equal to unity to better visually emphasize the shape 
and width of the recovered distribution. In the lower-left cor¬ 
ner of Fig. [3] we show the 2D joint marginalized likelihood 
distributions. In each panel, the cross denotes the location 
of our fiducial EoR parameters. For each, we construct both 
the la and 2cr contours (thick and thin lines, respectively), 
by computing a smoothed 2D histogram of the entire MCMC 
sample and estimating the likelihood contours which enclose 
68 and 95 per cent of the data sample, respectively. 

For all EoR parameters the recovered PDFs are cen- 
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Figure 3. The recovered constraints from 21CMMC on our three parameter EoR model parameters for a single (z = 9) 1000 h observation 
of the 21 cm PS obtained with HERA (red curve) and the SKA (blue curve). In the diagonal panels we provide the ID marginalized 
PDFs for each of our EoR model parameters (£o, R m fp and logio (7P| ed ) respectively) and we highlight our fiducial choice for each by the 
vertical dashed line. Additionally, we cast our ionizing efficiency, (fo. into a corresponding escape fraction, / esc , on the top axis (simply 
using the fiducial values in equation [2j. In the upper right panel we provide the ID PDF of the recovered IGM neutral fraction where 
the vertical dashed line corresponds to the neutral fraction of the mock 21 cm PS observation (igi = 0.71). Finally, in the lower left 
corner we provide the 1 (thick) and 2er (thin) 2D joint marginalized likelihood contours for our three EoR parameters (crosses denote 
their fiducial values, and the dot-dashed curves correspond to isocontours for Shi of 20, 40, 60 and 80 per cent from bottom to top). 


tred around their input fiducial parameters, highlighting 
the strong performance of 21CMMC at recovering our in¬ 
put model. This is despite the asymmetric distributions of 
the recovered PDFs, emphasizing the strength of the MCMC 
approach within 21CMMC. Comparing these two 21 cm ex¬ 
periments, we note the significant improvement achievable 
with the increased sensitivity of the SKA. This can easily 
be seen in the case of both £o and T \, where the widths 
of the ID PDFs, as well as the 2D joint likelihood contours, 
are noticeably tighter for SKA than for HERA. 

The improved constraints from the SKA relative to 
HERA can be best explained by referring to Fig. [2] We 
observe that while the sensitivity at large scales for both 
SKA and HERA are comparable, only the SKA is capable 
of probing the available small-scale information. Given the 
degenerate nature of the EoR parameters, by accessing this 
additional shape information tighter constraints are achiev¬ 
able. While it is not sufficient to break the degeneracies it 
does reduce their amplitude, substantially limiting the al¬ 
lowed EoR parameter space. It is important to note that 


with only a single redshift observation, the uncertainty is 
large and the likelihood is non-Gaussian distributed. Never¬ 
theless, our Bayesian MCMC framework captures parameter 
constraints and degeneracies, recovering the actual shape of 
the likelihood distribution (without the Gaussian assump¬ 
tions of the Fisher matrix formalism; Pober et al.|20141. 


In the £o-log 10 T(f i i. ed plane for HERA in lower-left panel 
of Fig. [3] we observe at the 95 per cent confidence level 
a multimodal distribution for £0 and T^ ed as highlighted 
by the lower streaks for high (0 and low Tfs cd . This region 
of parameter space is capable of reionizing the IGM earlier 
than our fiducial EoR model, due to a brighter population of 
lower mass galaxies (see Section [3X2}. This is highlighted 
by overlaying the xhi isocontours for the EoR parameter 
space, where it is clear that these two streaks reproduce dif¬ 
ferent IGM neutral fractions. Correspondingly, these models 
result in a smaller, secondary feature around Shi = 0.4 in 
the IGM neutral fraction PDF. Models in this region of the 
EoR parameter space have similar large-scale 21 cm power, 
but different small-scale structure which HERA cannot dif- 
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Instrument 


Parameter 


%I 

(single- 2 :) 

Co 

-^mfp 

l°gio(T F i f d ) 

z = 9 

HERA 331 

38.38±?2.43 

14.34+2;” 

4 tS7+°- 36 
-0.32 

0 7S+0- 10 
u - *°-0.21 

HERA 331 (with prior) 

38.371??;*® 

I3.88I4 39 

4 fi^+0.32 

4.00 — 0.35 

f) 7e;+0 09 
u - '°-0.13 

SKA 

09 1 1 +8.32 

OZ.±l_ 6< 47 

13.781*;*^ 

4.DU — 0.21 

0 71 

u - ‘ —0.09 


Table 2. The median recovered values (and associated 16th and 84th percentile errors) for our three EoR model parameters, £o> R m f P 
and T F ® ed and the associated IGM neutral fraction, Shi- We assume a total 1000 h integration time for a single epoch (z = 9) observation 
of the 21 cm PS and compare the recovered constraints for HERA, with and without a z = 7 IGM neutral fraction prior (Shi > 0.05), 
and for the SKA respectively. Our fiducial parameter set is (£o, i? m f p , logn)7^® ed ) = (30, 15 Mpc, 4.48) which results in an IGM neutral 
fraction of Shi = 0.71. 


ferentiate. The SKA, due to its higher sensitivity on small 
scales does not exhibit the same behaviour. Alternatively, 
these degeneracies can be broken with additional redshift 
observations which constrain the redshift evolution of the 
large-scale power (see below). 

In order to quantitatively assess the performance of 
21CMMC we choose to report the median value of each EoR 
parameter and the associated 16th and 84th percentiles. 
This choice follows on from the fact that the recovered 
ID PDFs do not need to be normally distributed within 
the Bayesian MCMC framework, as observed in the case of 
HER A in Fig. [3] In Table [2] we summarize the recovered 
EoR parameter constraints for each of the 21 cm experi¬ 
ments from a single epoch observation of the 21 cm PS. 

For HERA, we note that while both R m f p and T^ ed 
(also j:hi), are recovered to within 5 per cent of our fiducial 
parameters, in the case of £o we over predict the expected 
value by close to 30 per cent. This is not surprising given the 
mild positive asymmetry observed in the recovered ID PDF 
for Co- Since we choose to report the median recovered value 
for all EoR parameters, whenever we recover an asymmetric 
distribution our reported value can be shifted away from 
the input fiducial value. Despite the skewness, the 16th and 
84th percentiles should for the most part enclose our fiducial 
values. 

In the case of the SKA, we are able to recover all EoR 
parameters values at less than 10 per cent deviation from 
their fiducial values. While the errors on R m f p remain rela¬ 
tively consistent with those for HERA, we observe a reduc¬ 
tion of close to a factor of 2 (3) for the lower (upper) 1 <t 
bounds on Co and a 30 per cent reduction on the la errors 
on T F ? ed . We note a reduction also of approximately 25 per 
cent on the recovered 1 <t errors for the IGM neutral fraction. 

3.3 Single redshift observation with a neutral 
fraction prior 

In the previous section, we noted that the same large-scale 
power could occur at different IGM neutral fractions. This 
degeneracy can be broken with increased small-scale sensi¬ 
tivity (as in the case of SKA). Here, we investigate whether 
we can improve upon the available constraining power from 
a single redshift observation by instead including a prior on 
the IGM neutral fraction from observations at the tail end of 
the reionization epoch. For example, such constraints can be 
obtained from the spectra of high-z quasars such as ULAS 
<11120+0641 ( |Bolton et al. |[2011| |Mortlock et al.||2011| ) or 
from the observed drop in Lyman-a emission from high-z 


galaxies (e.g. Caruana et al.|2014 

Pentericci et al. 2014| Di- 

jkstra et al. 2014 Mesinger et al. 

2015). Here, we consider a 


purely illustrative conservative lower limit on the IGM neu¬ 
tral fraction at z = 7 of 5 per cent. We caution however that 
priors on different epochs can make our results more model 
dependent. For example, our fiducial EoR model does not 
allow for a redshift evolution in the reionization parameters. 
This is straightforward to account for, however, we postpone 
it to future work. 

In Fig. [4] we compare the recovered EoR parameter 
constraints for HERA from the same 1000 h observation of 
the 21 cm PS at z = 9, following the inclusion of the IGM 
neutral fraction prior at z = 7. Additionally, in Table [2] we 
again provide the recovered median and la errors for our 
EoR parameters. After the inclusion of the z = 7 neutral 
fraction prior we observe only a very marginal improvement 
in either the median recovered values or their 1 <t errors and 
this is equivalently shown by the slightly narrower widths of 
the recovered ID PDFs. 

However, in the case of the 2D joint likelihood distribu¬ 
tions we are able to observe an improvement in the allowed 
EoR parameter space. Specifically, focusing again on the £ 0 - 
log 10 T F ? ed plane, while the la errors are only marginally 
improved, the 2 <j contours are notably reduced. Most im¬ 
portantly, the IGM neutral fraction prior removes the sec¬ 
ond degenerate streak which marginally allowed a signifi¬ 
cantly different reionization history in the absence of the 
prior when observed with HERA. This can clearly be seen 
from the removal of the previously noted secondary bump 
at ihi = 0.4 in the IGM neutral fraction PDF. By imposing 
that the IGM must at least be 5 per cent neutral at z = 7, 
this rules out reionization scenarios which would otherwise 
predict an early end to the reionization epoch from a single 
redshift observation of the 21 cm PS. This highlights the 
importance of the additional leverage gained from probing 
the reionization history at improving the EoR reionization 
constraints, especially in the absence of strong telescope sen¬ 
sitivity at small scales. 

The level of improvement achievable from an IGM neu¬ 
tral fraction prior on a single epoch measurement of the 
21 cm PS depends on numerous factors. First, imposing a 
harder prior than our conservative choice would allow in¬ 
creased constraining power on the allowed EoR parameter 
space. Secondly, with decreasing telescope sensitivity, a neu¬ 
tral fraction prior becomes more important, as knowledge 
of the EoR history compensates for poor single epoch PS 
coverage. It is feasible that with improved foreground re¬ 
moval strategies, a marginal detection of the 21 cm PS could 
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<0 K mf p (Mpc) log 10 (T;ir d [K]) 


Figure 4. Similar to Fig. [ 3 ] except now we compare the impact of including an IGM neutral fraction prior for a single (z = 9) 1000 h 
observation with HERA. We consider an IGM neutral fraction prior at z = 7 of Xu , > 0.05 (a conservative choice motivated by recent 
QSO and LAE observations; green curves) and without a neutral fraction prior (red curve). The inclusion of a neutral fraction prior 
improves the constraints on the reionization parameters from a single z observation. 


be achieved with first generation instruments (Pober et al. 


2014). Therefore, in this scenario, a neutral fraction prior 


would be of vital importance for improving the available 
EoR constraining power. 

Finally, due to the flexibility of 21CMMC we emphasize 
that it is straightforward to additionally include direct priors 
on our model EoR parameters motivated by current and 
upcoming observations and theoretical advances. Given the 
observed degeneracies between the EoR parameters within 
this model, it is clear that including priors on our model 
parameters would notably improve the overall constraining 
power. 


3.4 Combining multiple redshift observations 

Here, we consider combining multiple epoch observations of 
the 21 cm PS, assuming a broader instantaneous bandwidth 
is available for the EoR observations. However, this comes at 
the cost of increased model dependence. For example, our 
three empirical EoR parameters are redshift-independent, 
whereas physically, we would expect each parameter to vary 
with redshift across the reionization epoch. In this sense, 
our astrophysical parameters can be thought of as ‘effec¬ 
tive’ EoR properties, averaging over any redshift evolution. 


In future, we will explicitly include a generalized redshift 
evolution for EoR parameters (e.g. Barkana 2009; Mesinger 
et al.|20T2 Zahn et al.||2012 l. 


To combine all the available information across the mul¬ 
tiple epochs and differentiate between the sampled EoR pa¬ 
rameters, we compute the total maximum likelihood fit of 
the EoR model, obtained by linearly combining the individ¬ 
ual x 2 statistics from each redshift. In principle, we could 
instead perform a weighted summation to estimate the to¬ 
tal maximum likelihood; however, in neglecting to do so our 
errors are more conservative. 

In Fig. [5] we compare the recovered ID marginalized 
PDFs and 2D marginalized likelihood contours for our three 
fiducial EoR parameters (Co, -R m f P and logio(T^r ))■ For 
this, we assume three independent 1000 h observations of the 
21 cm PS at z = 8, 9 and 10 for LOFAR (turquoise), HERA 
(red) and SKA (blue) each with a 1000 h integration time. In 
Fig ,[6]we provide the recovered ID PDFs of the IGM neutral 
fraction from each individual epoch. Note, for our fiducial 
EoR model, the corresponding IGM neutral fractions are 
Shi = 0.48, 0.71 and 0.83, respectively. Finally, in Table [3] 
we provide the recovered median, 16th and 84th percentiles 
for each EoR parameter and the corresponding IGM neutral 
fraction. 
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Instrument 

(multiz-z) 

Co 

Parameter 
Rmlp (Mpc) 

l°gio(T v !f d ) 

2 = 8 

x n i 

2 = 9 

2 = 10 

LOFAR 

HERA 331 

SKA 

45.21 ±?|;g5 
30.17±™° 
30.04t|;|| 

12-98 !£;" 
i5.i3±§;g? 
15.22±|;|i 

4.73l°;|| 

4.49±8;1§ 

4.48±°;" 

0 57+0-20 
u -°‘ -0.21 
f) 4Q+0-05 
U - 4y -0.05 

0 4Q+0-04 
U - 4y -0.04 

0 7fi+°-13 
u -' D —0.15 

n 71+0-04 

U - * -*--0.04 

0.7118:81 

o.8sl8:?8 

o.84i8:81 

0 S4+0-02 
U - y4 -0.02 

HERA 331 (w/o 25 per cent modelling uncertainty) 
SKA (w/o 25 per cent modelling uncertainty) 

28.78±?;|| 

30.31+1$ 

14.33±J;i| 

15.4711;" 

4 47+O.O6 
4>4 '-0.06 

4 4Q+0-04 
4 - 4 °— 0.04 

0.471q'q2 

o.48lg;81 

o.69i8:8i 

o.7ii8:81 

n ao+o.oi 
U.OZ_0.01 

n sQ+o.oi 
U.O4_0.01 


Table 3. Summary of the median recovered values (and associated 16th and 84th percentile errors) for our three EoR model parameters, 
Co? -R m f p and T^? ed and the associated IGM neutral fraction, xhi- We assume a total 1000 h integration time for each 21 cm PS 
observation at z = 8, 9 and 10 for LOFAR, HERA and the SKA, respectively. For both HERA and the SKA, we provide the recovered 
median values both including and excluding our 25 per cent modelling uncertainty (see Section |3.5| ). Our fiducial parameter set is (£o, 
R m f p , l°gio^Ji? ed ) = (30, 15 Mpc, 4.48) which results in an IGM neutral fraction of 5 hi — 0.48, 0.71, 0.83 at z = 8, 9 and 10 respectively. 
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Figure 5. The recovered constraints from 21CMMC on our reionization model parameters from combining three independent (z = 
8,9 and 10) 1000 h observations of the 21 cm PS. We compare the different telescope arrays, LOFAR (turquoise), HERA (red) and 
SKA (blue). Across the diagonal panels we provide the ID marginalized PDFs for each of our reionization parameters [£o, R m fp and 
logio(T^ ed ) respectively] and highlight our fiducial model parameter by a vertical dashed line. Additionally, we provide the 1 (thick) 
and 2cr (thin) 2D joint marginalized likelihood contours for our three reionization parameters in the lower left corner of the figure (crosses 
mark our fiducial reionization parameters). Including redshift information can significantly improve the constraints on each reionization 
parameter, and could even yield sufficient information to allow LOFAR to provide marginal constraints. 
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Figure 6. ID PDFs of the recovered IGM neutral fraction from our 1000 h observations of the 21 cm PS at z = 8 (left), z = 9 (centre) 
and z — 10 (right). Colours are the same as Fig.[5]and vertical dashed lines denote the neutral fraction of the mock 21 cm PS observations, 
Shi = 0.48, 0.71 and 0.83 respectively. 


From Fig. [5j we observe the ID PDFs and joint 
2D likelihood contours to be almost indistinguishable be¬ 
tween HERA and the SKA, with the SKA performing only 
marginally better. Furthermore, comparing with Fig. [3] the 
widths of both the ID PDFs and 2D likelihood contours 
have significantly reduced, drastically in the case of HERA. 
This emphasizes the strength of combining multiple epoch 
observations of the 21 cm PS (under the simple assumption 
they can be linearly combined). Quantitatively, we find lit¬ 
tle difference between the recovered constraints, with both 
instruments recovering median values of our EoR parame¬ 
ters to better than 1.5 per cent. Due to the improved con¬ 
straining power from sampling the reionization history the 
recovered PDFs are now normally distributed. Therefore, we 
can estimate the total ler fractional error for each parame¬ 
ter. For SKA (HERA), we find errors of 16.7 (22.0) on Co, 
17.8 (18.4) on R m f p and 2.4 (3.3) per cent on logio(T^® ed ). 
Additionally, across all three observed epochs, both HERA 
and SKA tightly recover the expected IGM neutral fraction 
of our fiducial EoR model (see Fig. |6|. 


Importantly, while for a single epoch observation, the 
SKA notably outperforms HERA, the improvements avail¬ 
able when probing the reionization history can make these 
two instruments comparable. By combining multiple epoch 
observations, these two experiments are able to pick up the 
same evolution in the 21 cm PS as the power rises and 
falls during reionization (see, e.g. Lidz et al. 2008 and our 
Section 3.1.21, compensating for their different sensitivity 
at small scales (for the case of these EoR models). There¬ 
fore, while the additional small-scale power can aid signifi¬ 
cantly in improving the constraining power from an individ¬ 
ual observation, it only marginally improves the constraining 
power across multiple epochs. 


For LOFAR, we find that by jointly combining mul¬ 
tiple epoch observations we are now capable of recover¬ 
ing reionization constraints. As expected though, the qual¬ 
ity of these constraints is drastically reduced compared 
to the second-generation experiments. Additionally, we re¬ 
cover broad marginalized ID PDFs, which encompass the 
majority of the EoR parameter space and their asymmet¬ 
ric nature complicates a quantitative comparison with the 
second-generation experiments. We find LOFAR can only 
marginally rule out (at ~ 2cr) either a rapid reionization 
(ihi < 0.2 at 2 = 8) or a delayed/extended reionization 
(ihi > 0.9 at a = 8). However, this is still encouraging for 


the first generation of 21 cm experiments given our conser¬ 
vative foreground and sensitivity assumptions. 


3.5 Impact of the EoR modelling uncertainty 


Thus far we have discussed the available constraints on our 
EoR parameters including a 25 per cent EoR modelling 
uncertainty, which dominates on large scales (see Fig. [2j|. 
Characterization of these uncertainties (e.g. by calibrating 
to more detailed RT simulations) will be essential to maxi¬ 
mizing the achievable scientific gains from upcoming second- 
generation interferometers. To emphasize this, in Fig. [7] we 
perform the same analysis as in the previous section, in¬ 
stead now removing the 25 per cent error for HERA and 
the SKA, corresponding to a perfect characterization of the 
EoR modelling uncertainties. We also list the correspond¬ 
ing constraints in Table [3] As expected, the resulting frac¬ 
tional errors are dramatically improved by a factor of 2-3 
to 8.1/5.5, 7.2/8.1 and 1.4/0.9 per cent for HERA/SKA (for 
Co, -Rmfp and log T^-® ed , respectively). 

The fact that the parameter constraints improve dra¬ 
matically and in a comparable manner for both HERA and 
the SKA shows that the redshift evolution of large-scale 
power dominates the constraining power for this EoR model. 
Unlike for narrow-band, single-z observations in Fig. [3] the 
additional leverage provided by small scales does not dra¬ 
matically improve the parameter constraints if the redshift 
evolution of the large-scale power can be constrained. We 
caution however that this conclusion is dependent on the 
EoR model. 

We also note that these values without the modelling 
error are in reasonable agreement with the rou ghly compa- 


Pober et al. 


(2014), sug- 


rablq | constraints recovered by 
gesting the Fisher matrix approach provides a reasonably 
accurate description of available constraints for this simple 
EoR parametrization. 


12 A like-to-like comparison is difficult for several reasons. First, 
these authors jointly combine four different epoch observations of 
the 21 cm PS compared to our three. Since instrument sensitivity 
increases with decreasing redshift, improved constraining power 
is available from their additional z = 7 observation (although it 
is not a significant improvement as three redshift measurements 
will have already strongly constrained the reionization history). 
Secondly, we consider only a 331 telescope antenna for HERA, 
unlike their 547 telescope design, while for the SKA there is the 
previously noted numerical error in their total sensitivities. 
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Figure 7. Analogous to Fig. [5] except here we show the impact of our assumed choice of a 25 per cent EoR modelling error, which 
dominates on large-scales (see Fig. [2j. Constraints from the SKA with the modelling error are shown in yellow, while those without 
any modelling error are shown in blue (SKA) and red (HERA). We see that under the optimistic assumption that we can perfectly 
characterize the EoR modelling uncertainties, the derived parameter constraints can be improved by a factor of a few (see Table [3|. 


4 TWO GALAXY POPULATIONS WITH 
MASS-DEPENDENT IONIZING 
EFFICIENCIES 

In the previous section, we highlighted the strength of 
21CMMC at providing constraints on our model parameters 
from a simple single ionizing source population. However, 
physically we expect galaxy formation to be more compli¬ 
cated. Therefore, in this section, we consider a more flexible 
EoR model which allows for both a halo mass-dependent 
ionizing efficiency and multiple ionizing source populations. 


4.1 Generalized five parameter EoR model 

Rather than assuming a step-function for the ionizing effi¬ 
ciency as in the previous section, here we allow it to vary 
with the host halo potential both for the bright, star-forming 
galaxies and the faint, feedback limited galaxies, 



r ( u ir y 

l ^Feed J 

T v i r 

^ rp Feed 
^ -^vir 


H 

( Y 

l fyiFeed J 

10 4 

K S) Tvir < T£r d 

(9) 


“0 

T v ir 

< 10 4 K. 



Here, £o, is simply the normalization of the ionization effi¬ 
ciency and retains the same definition as in Section ra 
When p = 0 and a —> oo, this expression is equivalent to 
the step function ionizing efficiency we considered earlier. A 
mass-dependent ionizing efficiency has previously been in¬ 
vestigated to describe the ionizing galaxy population; how¬ 
ever, these only considered a single population of galaxies 


described by a single power-law index (e.g. Furlanetto et al. 


2006a McQuirm et al. 2007 Paranjape & Cho udhury|2014 ). 
A sharp mass-scale separating feedback limited and non¬ 
limited galaxies is expected in some theoretical models (e.g. 
Finlator et al.|2011 Raicevic et al.|[2011 ). 

Our generalized EoR model contains a total of five pa- 
ramete rs, t he same three EoR parameters as described in 


Section 2.1 (£o, Rmfp and T/i® 8 ) and the two new power-law 


indices describing the mass-dependent ionizing efficiency, a 
and /?. For our fiducial mock observation, we choose to adopt 
a = 4/3 and f3 = —1/3, both within the allowed range, a, 
P £ [—4/3,10/3]. If instead we define equation § with re¬ 
spect to mass, these power-law indices become a = 2 and 
P = —1/2 (since M oc T^B). Our specific choice of a and 
P for our fiducial model is purely illustrative, and allows 
a strong mass-dependent scaling for the feedback limited 
galaxies, with a flatter slope contribution from the high-mass 
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Figure 8. A schematic picture of the z = 8 non-ionizing UV lumi¬ 
nosity function (blue curve) from our fiducial reionization model 
generalized to include two different source galaxy populations. 
First, above a critical temperature, T^? ed , haloes are sufficiently 
massive to be weakly affected by internal feedback process (e.g. 
supernovae driven winds), whose efficiency we describe by the 
power-law index /3. Between T^? ed and a cooling threshold tem¬ 
perature of T v i r = 10 4 K, haloes are sufficiently large to cool and 
form galaxies; however, their ionizing efficiency is affected by in¬ 
ternal feedback processes, described by a second power-law index 
a. Note that the model curve at T v [ r > T^? ed is obtained from 
a power-law fit to the halo mass as a function of UV magnitude 
recovered from abundance matching, which we can use to imply 


that for our fiducial choice of 0 = —1/3 we roughly obtain the 
scaling, f eBC oc M~°' 18 (see the text for further discussions). 


galaxies. In the local Universe, observations show that for 
low-mass galaxies the star formati on efficiency may increas e 
with host halo mass, /* oc M 2 / 3 (Kauffmann et al. 20031. 
On the other hand, our framework also allows small-mass 
galaxies to have even higher ionizing efficiencies than the 
high-mass galaxies. Such models (e.g. Alvarez et al. 2012) 
could be motivated by a much higher escape fraction in 
small galaxies (e.g. Wise et al. 2014). For our remaining 
EoR parameters, we adopt £o = 70, log 10 l(5)® ed = 4.7 and 
Rmfp = 15 Mpc. 

Again, we provide an illustrative example of our gener¬ 
alized EoR model by constructing a schematic picture of the 
expected z = 8 UV luminosity function in Fig. [8] following 
the same approach as in Sectio n |2.1| We convert our new 
critical feedback temperature, T v ? e into a UV magnitude 
(total halo mass) of Muv = —13.7 (M ~ 10 9 ' 11 Mq), and 
denote this and the atomic cooling threshold by the vertical 
dashed lines. For the massive (bright) galaxies, denoted by 
P = —1/3, we use the same fits to the UV luminosity func¬ 
tion in |Kuhlen fe Faucher-Giguere (20121. As before, this fit 
for these bright galaxies can be translated to / esc oc M -0 ' 18 
(see the discussion in Section 3.1.11. 


4.2 Constraints from combining multiple redshift 
observations 

Given the increased parameter degeneracies in this general¬ 
ized five parameter EoR model, we will focus only on mul¬ 


tiple epoch observations of the 21 cm PS. As in Section [3.4| 
we assume three independent 1000 h observations of the 21 
cm PS at 0 = 8, 9 and 10; however, only for HERA and 
the SKA. In Fig. [9] we show the recovered ID marginalized 
PDFs for each of the five EoR model parameters from HERA 
(red curves) and SKA (blue curves) and in Fig. 10 we show 
the corresponding 2D joint marginalized likelihood contours 
and the resulting ID PDFs of the IGM neutral fraction at 
each observed epoch. Finally, in Table [4] we report the re¬ 
covered median values and the corresponding 16th and 84th 
percentiles of our EoR model parameters. 

For the majority of these EoR parameters, we recover 
rather broad PDFs, noting mild to strong asymmetries em¬ 
phasizing the strength of the degeneracies. This is to be ex¬ 
pected with more model parameters. As an example, focus¬ 
ing on the a-log lo T(0® ed panel for a decreasing T^ ed (increas¬ 
ingly narrower mass range for the feedback limited galaxies), 
any strong suppression (high a) can mimic the same be¬ 
haviour as our fiducial model. This is due to the fact that 
for our assumed value of a = 4/3 the strong mass suppres¬ 
sion in the ionizing efficiency only allows low-mass galaxies 
closest to T^S cd to provide any meaningful contribution to 
reionization. 

Despite these degeneracies, we are still able to recover 
meaningful constraints on these EoR parameters. For exam¬ 
ple, the contribution of low-mass galaxies to reionization as 
well as a well-defined feedback scale are picked up in the 
marginalized ID T^ ed PDFs. Such fundamental issues in 
galaxy formation are unlikely to be addressed with direct 
observations, even with JWST. 

Quantitatively, both HERA and the SKA are able to 
recover median values for the five EoR parameters to within 
5 per cent of their fiducial values with a being the excep¬ 
tion (~ 20 per cent). The corresponding percentile errors 
for a are additionally rather broad, however, this was to be 
expected from the qualitative analysis above (i.e. a strong 
mass dependence of a is hard to constrain due to the nar¬ 
row range in mass of these galaxies which can contribute to 
reionization). Finally, from Fig. 10 we see that the PDFs 
of the IGM neutral fraction are tightly centred around the 
expected values from our fiducial model, and their quantita¬ 
tive errors are equivalent to our simple three parameter EoR 
model. This indicates that we have extracted all the avail¬ 
able information from sampling the reionization history for 
both our 21 cm experiments. In principle, we could improve 
these constraints by (i) adding measurements of the 21 cm 
PS in additional bands (e.g. at z = 7); (ii) applying priors on 
to our EoR model parameters; (iii) considering alternative 
statistics of the 21 cm signal in combination with the 21 cm 
PS; or (iv) reducing the intrinsic modelling uncertainty. 

It is important to note, that the discussion of the degen¬ 
eracies and constraints of these EoR parameters is specific to 
our chosen fiducial model for the mock observation. Due to 
the sharp break in the power-law indices describing the two 
populations of ionizing galaxies and our choice for T^ ed , 
the contribution to reionization from the high-mass (effi¬ 
cient star-forming) galaxy population can be relatively well 
constrained (as evident by the recovered PDF for /3). Such 
a model is motivated by theoretical predictions of a sharp 
feedback scale (e.g. Finlator et al. 2011 Raicevic et al.|20TT |. 
If however, we select a shallower break in the power-law in¬ 
dices, we effectively spread out the allowed range of masses 
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Instrument 
(multi- z) 

Co 

-^mfp (Mpc) 

Parameter 

a 

0 

logic (TFf d ) 

2 = 8 

®HI 

2 = 9 

2 = 10 

HERA 331 

SKA 

68.28 ±llf 6 

70.20±«;i 

14.99±l;|g 

14.87±i;H 

1 iq+1-28 

i - iU -0.72 

111 +1-^3 
1.11—0.70 

-o.3i±gi§ 

-0.33±g;g 

4-74±°;i 

4.73±°;l d 

0 28+ 0,04 
u - zo -0.03 

fl oq+0-03 
U.ZO — 003 

U.OO — 0.04 

u.OO_0.03 

0 74+0-04 

U - ' 4 -0.04 

0 74+0-03 
u -' 0.03 


Table 4. Summary of the median recovered values (and associated 16th and 84th percentile errors) for the five reionization parameters, 
Co, -R m fp, o, p and from the generalized EoR model (two ionizing galaxy populations) and the associated IGM neutral fraction, 

x H |. We assume a total 1000 h integration time for each of the three different epoch observations of the 21 cm PS at 2 = 8, 9 and 10 
for HERA and SKA. Our fiducial mock observation corresponds to (Co, R m fp, a, /!, logjQ'ZT® 6 ' 1 ) = (70, 15 Mpc, 4/3, —1/3, 4.7) which 
results in an IGM neutral fraction of ijji = 0.28, 0.55, 0.74 at z = 8, 9 and 10 respectively. 



Figure 9. The marginalized ID PDFs for each of the five reionization parameters (Co, -R m f p , <y - 0 out] T// e<l ) from our generalized EoR 
model containing two distinct ionizing galaxy populations for a combined 1000 h observation of the 21 cm PS across multiple redshifts 
(2 = 8, 9 and 10). We consider two separate telescope designs, HERA (red) and the SKA (blue), and the vertical lines denote our fiducial 
input parameters (Co, R m f p , a , 0 , T^? ed ) = (70, 15 Mpc, 4/3, -1/3, 4.7). 


of the ionizing galaxies (i.e. a shallower, positive a allows 
an increased contribution from the lower mass galaxies as 
they are not as efficiently suppressed). To illustrate this pes¬ 
simistic case, in Fig. [ill we consider (Co, -R m f P , a, /3, T^® ed ) 
= (30, 15 Mpc, 2/5, 1/15, 4.8). As expected, we now recover 
much broader PDFs. Most notably, with the shallower break 
in the power-law indices, it becomes increasingly difficult to 
constrain T^® ed due to the smoother transition between the 
ionizing efficiency of the two galaxy populations. This is un¬ 
derstandable, as in this model the two galaxy populations 
resemble a single population. 


5 CONCLUSION 

The EoR is rich in astrophysical information, probing the 
nature of the first stars and galaxies, as well as their impact 
on the IGM. However, our current understanding of EoR 
astrophysics is poorly understood. Our understanding will 
rapidly improve in the near future with the detection of the 
cosmic 21 cm signal from a broad range of dedicated 21 cm 


experiments. However, we are faced with the fundamental 
question, what exactly can we learn from these observations? 

To this end, we developed 21CMMCp^| a new, massively 
parallel MCMC analysis tool specifically designed to recover 
constraints on any astrophysical parameters during the EoR 
from observations of the cosmic 21 cm signal. 21CMMC uti¬ 
lizes existing astrophysical simulation codes and parameter 
estimation techniques to ensure its efficiency and flexibil¬ 
ity. It combines the astrophysical parameter sampler COS- 
MOHAMMER (Akeret et a.1.12013), calling the seminumerical 


21 cm simulation code 21CMFAST (Mesinger & Furlanetto 


2007 Mesinger et al. 20111 at each point in the sampled 


EoR astrophysical parameter space. This allows sampling 
over 10 5 individual realizations of the cosmic 21 cm signal 
to recover robust EoR parameter constraints, without com¬ 
mon assumptions of Gaussian errors. 

Throughout this work, we have emphasized several key 
strengths of 21CMMC, outlining its applicability to a broad 


13 http://homepage.sns.it/mesinger/21CMMC.html 
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Figure 10. The 2D joint marginalized likelihood contours for our generalized five parameter EoR model (£o, R m f p , ct, 3 and TT® ed ) for 
a combined 1000 h observation of the 21 cm PS across multiple redshifts (z = 8.9 and 10). We consider two separate telescope designs, 
HERA (red) and SKA (blue), crosses denote the fiducial model parameters and the la and 2cr contours are shown as the thick and thin 
contours respectively. Additionally, in the top-right corner we provide the one-dimensional PDFs of the recovered IGM neutral fraction 
for each individual redshift observation, with the vertical dashed line denoting the neutral fraction of the mock 21 cm PS observations, 
Shi = 0.28, 0.55 and 0.74 respectively. 


range of studies of the 21 cm signal during the EoR. In 
summary, we 


i) showed by utilizing a Bayesian MCMC framework, we 
can efficiently and accurately recover robust astrophysical 
constraints by recovering the marginalized ID PDFs, irre¬ 
spective of any model degeneracies. This removes the fun¬ 
damental limitations and assumptions inherent in previ¬ 
ous EoR astrophysical parameter studies such as sampling 
fixed grids (e.g. Choudhury fc Ferrara||2005 Barkana||2009 


Zahn et al. ||2012 ), Fisher Matrix ap 


Mesinger et al. 2012 


proaches (Pober et al. 


2014 1 or simplistic analytic modelling 


of the EoR ( Harker et al.||2(5T2 Morandi & Barkana 2012 


Patil et al. 20141; 
ii) highlighted that it is generalizable to any theoretical 
model of the EoR. We showcased two such models: (i) a 
popular, three parameter EoR model driven by a single pop¬ 
ulation of star-forming galaxies characterized by a constant 
(mass-independent) ionizing efficiency; and (ii) a generalized 
five parameter EoR model containing two ionizing galaxy 


populations each characterized by a unique mass-dependent 
power-law ionizing efficiency; 

in) highlighted the available constraining power from both 
single redshift and combining multiple redshift observations 
of the ionization 21 cm PS for three current and proposed 
21 cm experiments, LOFAR, HERA and the SKA; 

iv ) recovered for LOFAR/HERA/SKA fractional errors of 
45.3/22.0/16.7, 33.5/18.4/17.8 and 6.3/3.3/2.4 per cent on 
the ionizing efficiency (escape fraction), mean free path of 
ionizing photons and the minimum halo mass hosting star- 
formation from combining three independent lOOOhr obser¬ 
vations of the 21 cm PS at z = 8, 9 and 10 assuming no 
priors and a conservative noise estimate; 

v) showed that by removing the EoR modelling uncer¬ 
tainty from our fiducial analysis, these constraints can be 
improved by up to a factor of a ~few. This result moti¬ 
vates further work on characterizing EoR modelling errors 
in order to maximize science returns from second-generation 
instruments; 

vi) quantified the improvements in constraining power on 
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Figure 11. Same as Fig. [9] except now we consider an EoR model with a weaker break in the power-law indices. For this model, our 
EoR parameters are (£o> /i' ln f p , a, 0, log TP; ed ) = (30, 15 Mpc, 2/5, 1/15, 4.8). As expected, reducing the amplitude of each power-law 
index and weakening the relative break between them, reduces the overall constraining power on all EoR model parameters. 


the EoR parameters from accessing the additional informa¬ 
tion from the cosmic reionization history through the ap¬ 
plication of physically motivated priors on the IGM neutral 
fraction at the tail-end of reionization; 
vii) emphasized that 21CMMC additionally allows for the 
inclusion of individual priors on the EoR model parameters 
motivated by current and upcoming observations and theo¬ 
retical advances. Furthermore, while we have focused solely 
on the 21 cm PS, it is easily adaptable to consider any sta¬ 
tistical measure of the cosmic 21 cm signal. 

Finally, due its applicability to a broad range of EoR 
studies, 21CMMC could be an important analysis tool for 
the 21 cm EoR community. For example, it can be used 
to quantify how foregrounds and other contaminants can 
hinder the recovery of EoR astrophysics. Moreover, it can 
serve to guide designs and observing strategies for current 
and future 21 cm experiments to maximize their scientific 
returns. 
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